//1D PIC code to solve plasma two-stream instability problem

#include "mine.h"

double L;
int N,J;

int main() 
{
    // Parameters
    L;   // Domain of solution 0<=x<=L(in Debye lengths)
    N;   // Number of electrons
    J;   // Number of grid points
    double vb;   // Beam velocity
    double dt;   // Time-step (in inverse plasma frequencies)
    double tmax;   // Simulation run from t=0. to t=tmax

    // Get parameters
    FILE* input;
    input = fopen("input","r");
    if ( input == NULL ) {
	printf ( "Error - invalid opening file input\n" );
	exit (1);
    }
    fscanf(input, "%d,%lf,%lf,%d,%lf,%lf",&N,&vb,&L,&J,&dt,&tmax);
    int skip = int (tmax / dt) / 10;
    if ( (N < 1) || (J < 2) || (L<=0.0) || (vb <=0.0) || (dt <= 0.0) || (tmax <= 0.0) || (skip < 1) ) {
	printf ( "Error - invalid input parameters\n" );
	exit(2);
    } 
    else {
	printf ( "The number of electrons: %d\n", N );
	printf ( "The velocity of beam: %f\n", vb );
	printf ( "The length of simulation domain: %f\n", L );
	printf ( "The number of grid points: %d\n", J );
	printf ( "The time step: %f\n", dt );
	printf ( "The period of simulation time: %f\n", tmax );
    }

    //Set names of output files
    char* phase[11];
    char* data[11];
    phase[0] = "phase0.out"; phase[1] = "phase1.out"; phase[2] = "phase2.out"; phase[3] = "phase3.out"; phase[4] = "phase4.out"; phase[5] = "phase5.out"; phase[6] = "phase6.out"; phase[7] = "phase7.out"; phase[8] = "phase8.out"; phase[9] = "phase9.out"; phase[10] = "phase10.out"; 
    data[0] = "data0.out"; data[1] = "data1.out"; data[2] = "data2.out"; data[3] = "data3.out"; data[4] = "data4.out"; data[5] = "data5.out"; data[6] = "data6.out"; data[7] = "data7.out"; data[8] = "data8.out"; data[9] = "data9.out"; data[10] = "data10.out";

    // Initialize solution
    double t = 0;
    int seed = time (NULL); srand(seed);
    Array<double,1> r(N), v(N);
    Array<double,1> ne(J), rho(J), phi(J), E(J);
    for (int i = 0; i < N; i++) {
	r(i) = L * double (rand()) / double (RAND_MAX);
	v(i) = distribution (vb);
    }
    Output2 ( phase[0], t, r, v );

    // Evolve solution
    Array<double,1> y(2*N);
    for (int k=1; k <= 10; k++) {
	for (int kk=0; kk < skip; kk++) {
	    
	    // Calculate electron number density
	    Density (r, ne);
	    
	    // Solve Poisson's equation
	    double n0 = double (N) / L;
	    for (int j = 0; j < J; j++) 
		rho(j) = ne(j) / n0 - 1.0;
	    double kappa = 2.0 * M_PI / L;
	    Poisson1D (phi, rho, kappa);
		
	    // Calculate electric field
	    Electric (phi, E);

	    Load(r, v, y);
  
	    // Take time-step and develop the time
	    rk4_fixed (t, y, E, rhs_eval, dt);
	    
	    // Make sure all coordinates in range 0 to L
	    for (int i = 0; i < N ; i++) {
		if (y(i) < 0.0) y(i) += L;
		if (y(i) > L) y(i) -= L ;
	    } 
	    
	    // UnLoad the phase data seperatively
	    UnLoad (y, r, v);
	}
	printf("t = %11.4e\n", t);
	
	// Output data
	Output2(phase[k], t, r, v);
	Output3(data[k], t, ne, rho, E);
    }
	
    return 0;
}

